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Abstract 

A thin shell finite element approach based on Loop's subdivision surfaces is proposed, capable 
of dealing with large deformations and anisotropic growth. To this end, the Kir chhoff- Love theory 
of thin shells is derived and extended to allow for arbitrary in-plane growth. The simplicity and 
computational efficiency of the subdivision thin shell elements is outstanding, which is demonstrated 
on a few standard loading benchmarks. With this powerful tool at hand, we demonstrate the broad 
range of possible applications by numerical solution of several growth scenarios, ranging from the 
uniform growth of a sphere, to boundary instabilities induced by large anisotropic growth. Finally, it 
is shown that the problem of a slowly and uniformly growing sheet confined in a fixed hollow sphere 
is equivalent to the inverse process where a sheet of fixed size is slowly crumpled in a shrinking hollow 
sphere in the frictionless, quasi-static, elastic limit. 

1 Introduction 

Thin shells are vital parts of innumerable problems in structural engineering and material science. It is 
impossible to imagine many of today's technological achievements without a deep understanding of their 
large deformation response to external loads, spatial constraints, nonlinear material behavior, self-contact 
and other multifarious interactions. For solving the complex interplay of these effects, one resorts to 
numerical methods, like for the wrinkling of metal sheets in a vehicle crash (e.g., [1-4]), or the crumpling 
of paper [5-11], just to name a few. The finite element method (FEM) has proven to be among the 
most flexible and efficient tools for a large number of such problems, in particular where complicated 
geometries, strong material nonlinearities, or anisotropy come into play. 

The nontrivial nature of the various large strain modes of shells is due to their thinness, enabling 
a large variety of complex three-dimensional deformations induced by the aforementioned external or 
intrinsic constraints. In the KirchhofT-Love theory [12], which has been well understood and widely 
applied for decades, the thinness of the shell manifests itself in the assumption that material lines that 
are straight and perpendicular to the middle surface retain these properties and their length during shell 
deformations. To account for out-of-plane bending stiffness, the resulting total elastic energy formulation 
integrates the mean and Gaussian curvatures over the shell's middle surface. In the context of a finite 
element treatment of structural shell analysis, this inevitably calls for shape functions with continuous 
first derivatives (C 1 continuity), i.e., functions that belong to the Sobolev space H 2 . This requirement has 
posed a tough challenge in the history of shell finite elements, giving rise to a vast amount of literature 
on ways of coping with the disadvantageous consequences such as locking phenomena, resulting from 
insufficiently meeting this requirement. For a brief overview, see e.g. Ref. [13]. Many of the developed 
workaround methods introduce interpolation coefficients for higher derivatives of the displacement field, 
leading to a significant increase in the number of degrees of freedom (DOFs) to compute. Parisch [14] 
has proposed quadrilateral shell elements with only displacement variables, which however still require 
special measures against severe locking. An entire class of quadrilateral C 1 elements is obtained from 
tensor products of Hermite polynomials, see e.g. Ref. [15]. However, while simple to construct, they are 
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limited to regular rectangular meshes. The family of Hsieh-Clough- Tocher triangles have shown success 
in plate bending and other biharmonic problems [16], but are tedious to set up and add many additional 
DOFs. The recently developed subdivision surface shape functions [17,18], which strictly belong to the 
class C 1 everywhere in the domain, has successfully vanquished all these problems, and is hence gaining 
increasing interest in the FE community. 

External loads and constraints are not the only causes for nonlinear shell deformations. Thinking of 
nature's soft tissues such as leaves, flower petals, cell membranes, insect wings, etc., it becomes evident 
that large deformations of shells are often on account of growth or shrinkage (e.g., [19,20]). We use the 
term growth to represent both these mutually inverse processes in the following. Growth often leads to 
the inevitable development of residual stresses (e.g., [21,22]), causing a shell to deform in an attempt 
to minimize them. The study of growing thin shells is, however, not limited to bioengineering. There 
is also a large potential in the blossoming fields of bionics or material engineering, where a variety of 
smart or self-actuating materials is designed (e.g., [23,24]). Predicting the large deformation response 
of such growing thin sheets, be it for whatever physical cause or technical purpose, calls for an efficient 
and robust numerical tool inherently featuring the capability to grow according to arbitrary prescribed 
anisotropic growth fields. Not many common numerical discretization techniques are fit for this kind of 
anisotropy. Discrete elements, or beam networks, for instance, are not well suited, as their only degree of 
freedom capable of accounting for in-plane growth is the edge length connecting the mesh vertices, which 
may not be aligned with the desired growth direction. To the finite element method, on the other hand, 
material anisotropy poses no problem, since the mesh faces are numerically integrated over, regardless of 
their orientation. 

The purpose of this article is to present easy-to-implement thin shell finite elements that are capable 
of anisotropic in-plane growth, while profiting from the superior efficiency and strong robustness of 
the subdivision surface paradigm at the same time. The adopted continuum model is based on the 
geometrically nonlinear Kirchhoff-Love shell theory, which we extend by a large-strain continuum growth 
model that supports anisotropic growth fields. The basic idea is to assume that body deformations can 
be due to both a change of mass or volume and an elastic response [22,25]. The thin shell growth model 
employed here is based on the volumetric growth assumptions of Rodriguez et al. [26] , which have later 
been put on rigorous foundation [27,28]. 

The paper is organized as follows: The next section extends the Kirchhoff-Love theory of thin shells for 
large deformations by Rodriguez' growth ansatz, for the continuum shell and in discretized form suitable 
to finite element analysis. A short overview of subdivision surface interpolation for finite elements and 
boundary conditions is given in the subsequent section 3, followed by a description of some implementation 
details in section 4. Standard load cases to verify and benchmark our shell elements with and without 
growth are presented in the last section, where we also demonstrate the huge potential in applicability of 
growing thin shells to several problems in material science and engineering. 

2 The Kirchhoff-Love Theory with in-Plane Growth 

In the following section, the Kirchhoff-Love theory is briefly derived following the common stress-resultant 
formulation, meaning that the stresses are integrated analytically through the thickness, so that one is 
left with a resultant stress on the middle-surface. In the course, the theory is amended by anisotropic 
in-plane growth for large strains. 

2.1 Kinematics of Deformation 

Let Greek indices a,/3,7, S take the values 1 and 2, and Latin indices z,j take values from 1 to 3. Unless 
otherwise mentioned, the Einstein summation convention applies to repeated indices, and lower (upper) 
indices denote the covariant (contravariant) components. 

Let Q C E 3 be the geometry of the stress-free undeformed middle surface of a shell with small thickness 
/i, embedded in three-dimensional Euclidean space. Under the action of external loads or growth, the 
shell deforms into a new configuration characterized by the middle surface OcE 3 . Let {fl 1 ,^ 2 ,^ 3 } be 
a curvilinear coordinate system, and x(^ 1 ,6 >2 ) and x^ 1 ,^ 2 ) be parametrizations of the reference middle 
surface ft and deformed middle surface respectively, see Fig. 1. The positions r and r of material 
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points in the reference and deformed shell may be parametrized as 

r(9 1 ,9 2 ,6 3 )=x(6 1 ,6 2 ) + 6 3 H 3 (6 1 ,6 2 ), 6 3 G [—h/2, h/2], 

r(e 1 ,e 2 ,e 3 ) = x(e\e 2 ) + e 3 a 3 (e\e 2 ), e 3 e \-h/2,h/2). 

The tangent space of the middle surface is spanned by the vectors 



<9x 

00" 



and 8L a (6\6 2 ) = 



<9x 

36" 



(1) 
(2) 



(3) 



and by virtue of the Kirchhoff kinematic assumption, the material orientation in the thickness direction 
of the shell is determined by the shell directors 



a 3 



a x x a 2 
|ai x a 2 | 



and 



a 3 



ai x a 2 
|ai x a 2 | 



(4) 




Figure 1: Reference and deformed configurations of the shell middle surface with parameterizations 
x(6> 1 ,6> 2 ) and x((9 1 ,6> 2 ), respectively. 



The covariant components of the surface metric tensors, or first fundamental forms, follow as 

= a a • a^, a a (3 = & a • a^, (5) 

those of the shape tensors, or second fundamental forms, as 

K a p = -a 3 , a • a^ = a 3 • a^, « a/9 = -a 3 ^ • a^ = a 3 • a a ^, (6) 

and the infinitesimal area element can be expressed as dft = |a~i x a 2 | d6 1 d6 2 . The covariant basis vectors 
for a generic point within the shell are given by 



dr 

dO* 
dr 

dO* 



S a = = a « +6> a 3,«, 

Q3 r 



dr 
dr 



S3 ^ 3 a 3, 



and the covariant components of the corresponding metric tensors g, g are 

9ij = Si ' §j? 9ij ~ Si ' Sj- 



(7) 
(8) 

(9) 



Owing to the Kirchhoff constraints, Eq. (4), we require that g 33 = 1 and g a s = g^ a — 0. The deformation 
gradient in curvilinear coordinates reads 



F := V F r = 



dr 



1 S = Si ® s 



and maps between deformed and reference metric, 



(10) 



(11) 
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In the growth model proposed by Rodriguez et al. [26], the geometric deformation gradient is multi- 
plicatively decomposed into a growth tensor G and a purely elastic response A, that ensures compatibility 
and continuity of the body, according to 

F = AG. (12) 

We further assume the Green-Lagrange strain tensor E to depend solely on the elastic response, and G 
to be independent of the deformed configuration, such that 

E := \{g -g) = \(A T gA -g) = ^(G^G" 1 -g) = \{g- g), (13) 

where g := G -T gG -1 is the growth-modified metric. For the following, we restrict G to anisotropic 
in-plane growth, i.e., 

'[G afi ] 



O 1 1 



(14) 



Consequently, §33 = 1 and g a s = §3 a = 0, and the remaining components can be expanded in terms of 
the thickness parameter 3 to 

g a p = {[G a p]- T [a af} ] [Gcp]- 1 )^ - 20 3 ([G a/3 ]" T [ Kafl ] [G^]" 1 )^ + O((0 3 ) 2 ) (15) 
= a aP - 28 3 k aP + O((0 3 ) 2 ), (16) 

Analogous to the derivation of the Kirchhoff shell without growth, we neglect terms O((0 3 ) 2 ) in the 
following, making the theory a first order shell theory which is valid for small thicknesses h. The non- 
zero components of the Green-Lagrange strain (13) then follow as 

« a a p + 6 3 p a p (17) 

where a and j3 are the growth- modified membrane and bending strain tensors, respectively, 

OLot^ = ^iflotp ~ Q>otp) (18) 

fia(3 = K a p - R a {3. (19) 

The above expressions are formally identical to the strains of a Kirchhoff shell without growth. Note 
also that the four in-plane growth parameters G a £, which in general may be functions of various external 
or internal variables such as time, space, stress, etc. [29,30], are expressed with respect to the reference 
tangent basis {ai,a 2 } in the above formalism. In practice, it is thus necessary to perform a change of 
basis when they are to be given with respect to a specific coordinate system, such as the Cartesian, for 
instance. 



2.2 Constitutive Model 

Assuming that the shell obeys the St. Venant-Kirchhoff law, the connection between its geometrical 
configuration and material properties is provided by the Koiter energy density functional [31,32] 

W = ±KH a W& a p& 7S + \DH^ 5 ^P l5 (20) 

with membrane stiffness K and bending rigidity D, given by 

Yh Yh 3 
1 — v l 12(1 — v l ) 

and with the elasticity tensor 

H a ^ s = va^a? 5 + ^^(a a ^ 5 + a aS a^). (22) 

Y is the Young's modulus and v the Poisson ratio. The growth-modified resultant membrane stresses ft 
and resultant bending stresses rh follow by the principle of work conjugacy: 

= = KH a W&~6, rh a P = = DH a ^ 5 p l5 . (23) 

Oa a(3 d(3 a(3 

Notice that the energy density in the above form coincides with the non-Euclidean plate approach 
derived in Ref. [33]. 
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2.3 Variational Formulation 



The total potential energy <1> of the Koiter shell with total Lagrangian displacement of the middle surface, 
u = x — x, is obtained by adding the internal elastic energy <E> mt to the contribution <£ ext from external 
loads q per unit surface area and traction N per unit edge length, yielding 

*[u] = $ int [u]+$ ext [u], (24) 

$ int [u] = (_Wdti, (25) 
Jq 

$ ext [u] = - /qudft- [_N-uds. (26) 
Jo. Joq 

Our aim is to find the minimum of Eq. (24) for prescribed growth tensors, which is equivalent to 
finding a displacement field u satisfying the variational problem 

= 5<S>[u] = 5& nt [u] + £$ ext [u], (27) 

5& nt [u] = J_ (n a(3 Sa af3 + m a(3 5f} a p) dfi, (28) 

£$ ext [u] = - [qSudH- [ N • Su ds. (29) 
Jq Joq 

We further augment Eq. (27) with the usual inertial term to capture the dynamics of the system. The 
variational statement thus becomes 

= 5$ [u] + _ hpu • Su dO, (30) 
Jq 

where p is the mass density of the shell. The variation of the membrane and bending strains are easily 
calculated if the growth tensor is assumed not to depend on the displacement field: 

[S&ap] = [G a p]- T [Sa a p] [G af} ]-\ [5p al3 ] = [G a p]~ T [Sn a p] [G^]" 1 . (31) 

5a a f3 and SK a p are the usual first strain variations for shells without growth, see e.g. Ref. [34,35]. Since 
we will integrate the Newtonian equations of motion induced by Eq. (30) explicitly in time (cf. Section 
4.1), no second variations are needed at this point. We note, however, that the derivation of the second 
variations straightforwardly follows the usual formalism (e.g., [35,36]) without complications for our 
growth-modified strains. 

2.4 Finite Element Discretization 

The finite element discretization is done in the usual way: The minimization problem is replaced by an 
approximate minimization problem over a finite subspace Vh C V := i^ 2 (^,E 3 ) of admissible displace- 
ments: 

inf $[u] — > min $[ U/l ]. (32) 
uGV u h ev h 

Vh is spanned by a finite set of basis functions {Ni(6 1 ,# 2 ), / = 1, ...,A/" n } with local support, where J\f n 
is the number of mesh nodes. The displacement field is then written as a linear combination of the trial 
space basis functions: 

N N 
u h (0\0 2 ) = Y,"iNi(0\0 2 ), Su h (9\9 2 ) = '£5u I N I (0\0 2 ). (33) 
i=i i=i 

Substituting the above interpolation into the weak form, Eq. (30), and using the arbitrariness of the trial 
field, the variational statement is recast into an algebraic minimization problem for the nodal displace- 
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merits u/: 



= f} nt - ff xt 

f int = _ f 

Jn 

f'ext 
/ 

Mjj 



J 



/q7V 7 dn+ [_NN I ds J 

[hpNiNjdQ. 
Jn 



(34) 

(35) 

(36) 
(37) 



As usual in finite element analysis, the integrals (35-37) are evaluated numerically using a quadra- 
ture rule with J\f p integration points {q p = (0p,0p),p = 1,...,A/" P } and corresponding weights {w pi p = 
1, A/" p }, taking advantage of the local support of the shape functions. A single element's contribution 
to the generalized internal force (35), for instance, becomes 



r int 



<9uj 



ai x a 2 



(38) 



where [*](0i,02) denotes evaluation of the integrand at the quadrature point q p mapped onto element e. 
Analogous to Eq. (31), the growth- modified strain derivatives are straightforward to calculate since the 
growth tensor is independent of the deformed configuration: 



da 



a(3 



dur 



[G Q 



da 



a(3 



dur 



[G Q 



a/3 



dur 



[Ga£ 



dfr 



dur 



[G a pY 



(39) 



The energy functional (20) and the generalized internal force (35) contain second derivatives of the 
displacement field u. For boundedness of these integrals, the shape functions Nj therefore need C 1 
continuity. Loop subdivision surfaces, being C 1 -continuous everywhere and even C 2 -continuous except 
at a finite set of extraordinary points, fully meet this requirement. 

A noteworthy difference of the above finite element description of growth to recent approaches via the 
introduction of prescribed non-Euclidean target metrics [37-40] is that the numerical implementation of 
the present model includes the change of reference curvature when the surface grows, which is missing in 
the tethered mass-spring model of the target metric approach [37,40]. 



3 Subdivision Surfaces 

Subdivision surfaces were developed simultaneously by Catmull and Clark [41] in the context of computer 
graphics in 1978 as a method of representing smooth surfaces by a coarse polygonal mesh, termed the 
control mesh. Cirak et al. [17, 18] have ported them to the finite element method. 

3.1 Shape Functions 

The methodology is based on Stam's eigenanalysis [42] of Loop's recursive refinement rule [43] for tri- 
angulated surfaces with arbitrary topology, which gave access to a set of 12 quartic box splines, that 
exactly interpolate the infinitely refined surface, called limit surface, at all points except a finite set. 
In fundamental difference to traditional finite elements, subdivision surfaces gain C 1 continuity, or H 2 
integrability, required for a finite deformation energy of a shell, at the expense of a larger local support 
of the basis functions. Instead of only the 1-ring consisting of directly adjacent neighbor elements, each 
element's support spans also the 2-ring of next-neighboring elements. While this is neither a mathe- 
matical problem nor computationally disadvantageous, it does require new concepts for, e.g., boundary 
conditions, as explained in the next subsection. Details on the Loop subdivision shape functions and 
their application to interpolating the limit surface on arbitrarily triangulated meshes can be found in 
Refs. [17,42]. 

Subdivision surface elements are much more efficient than conventional C° shell elements, as will be 
demonstrated toward the end of this paper. Not only are they completely free of membrane and shear 
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locking phenomena, due to the fact that they fully meet the continuity requirements imposed by the 
continuum theory, negating the need to incorporate numerical workarounds for such problems. They also 
go without rotational degrees of freedom, i.e., each mesh node carries three displacement variables only. 
Moreover, a single Gauss point per element is sufficient for high accuracy in many cases [17,44]. 

3.2 Boundary Conditions 

The treatment of domain boundaries is in general non-trivial for subdivision surface interpolation because 
elements along a shell boundary lack a complete 1-ring and hence per se cannot be interpolated like 
elements in the interior of the domain. Three principal approaches have been proposed to solve this 
problem. 

1. As pointed out by Cirak et al. and Biermann et al. [17,45,46], a special subdivision rule may 
be applied to boundary elements, corresponding to one-dimensional subdivision on the boundary 
curve. 

2. Alternatively, Schweitzer [47] suggested appending an additional row of "ghost" vertices and el- 
ements to the control mesh along its boundaries. If these ghost vertices are positioned and the 
displacement field is projected onto them component-wise according to Cirak et al. [17], see Fig. 2, 
the usual boundary conditions can easily be imposed. We will refer to this type of boundary as the 
Schweitzer- Cirak type in the following. However, as demonstrated by Green [36,48], such boundary 
constraints are overly restrictive and lead to a reduction of the convergence with the number of ele- 
ments from order two to one. In practice, this implies that rather fine meshes are required for high 
accuracy, undermining the otherwise high computational efficiency of subdivision finite elements. 

3. The third approach was proposed by Green [36,48] as a remedy to the limitations of the second. 
Instead of drastically constraining the ghost displacements according to Fig. 2, only the minimum 
set of necessary conditions is imposed directly on the limit surface. The resulting linear constraint 
equations can be solved using the penalty method, Lagrange multipliers, or any other solving 
technique suitable for constrained minimization. 



BC type imposed displacements 

free u i =112 + 113 — ui 

pinned u 2 = u 3 = 0, u 4 = — ui 

clamped ui = 112 = 113 = U] = 



Figure 2: Boundary conditions of the Schweitzer- Cirak type. 

In summary, Schweitzer- Cirak boundaries are the easiest to implement, but should be avoided in cases 
where physical accuracy is crucial. The two alternatives are significantly more complex to implement in 
general, with the exception of free Green boundaries, where the ghost nodes are simply left unconstrained. 
Such ghost nodes are, however, unsuited for static analysis due to underdetermination of the system. Most 
boundaries in the numerical examples presented in Section 5 are free, and since we solve the dynamic 
problem, Green's convenient method is applied there. We use Schweitzer-Cirak boundaries in all other 
situations for simplicity. 

4 Numerical Implementation 

Thin shells are commonly formulated and implemented using Voigt's vector notation for convenience and 
efficiency, exploiting the symmetry of the involved tensors. Readers interested in the details are referred 
to the numerous publications describing the technicalities, such as Refs. [17,35,49]. The same principles 
apply straightforwardly to the growth- modified strains and stress resultants. In this section, we focus on 
the integration of dynamics and shell contacts. 
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4.1 Integration in Time 

We solve the hyperbolic equilibrium equations 

Mii(t) + Cu(t) + f int (u(t)) = f ext (u(t),t) . (40) 

Here, u(t) is the displacement vector containing the nodal interpolants u/ at time t, M is the mass matrix 
assembled from element contributions to Eq. (37), and C is a viscous damping matrix for equilibration. 
f mt and f ext account for internal out-of-equilibrium forces, Eq. (38), and external loads and contact forces, 
respectively. Newmark's family of integration methods [50] is widely used in structural dynamics [51,52]. 
Let u t w u(t), v t & u(t) and a £ « ii(t) be the discretized approximations of the displacement vector and 
its time derivatives. For fixed Newmark scheme parameters j3 and 7, they are integrated according to 

(At) 2 

Ut+At = u £ + At v £ + ((1 - 2/3)a £ + 2/3 a £+A£ ) , (41) 

Vt+At = v £ + At((l - 7 )a t + 7a t +At). (42) 

At denotes the finite time step. We apply the unconditionally stable constant- average acceleration method, 
that is obtained by setting f3 = 1/4 and 7 = 1/2, in form of a predictor-corrector scheme with lumped 
masses and subcritical lumped damping. For adaptive time step control, the a posteriori local error 
estimator by Zienkiewicz and Xie [51,53] is employed. 

4.2 Shell-Shell Contact 

For problems where the shell may be in contact with itself, we use a hierarchical spatial decomposition 
to find contact points efficiently, similarly to Ref. [54]. On the coarsest level, elements are placed into an 
array of cubic axis-aligned cells. Each cell holds a reference to all elements whose axis-aligned boundary 
boxes (AABBs), extended by half the shell thickness in each direction, overlap with it. Each element 
pair sharing at least one cell enters the mid-level check where the elements' extended AABBs are tested 
for overlap. If they do, the distance between the two elements and the corresponding closes points are 
computed on the finest level of the hierarchy. 

Finding the closest points between the subdivision limit surfaces of two triangles is a four-dimensional 
nonlinear optimization problem with linear inequality constraints. As such, it may be solved with, e.g., 
the projected gradient method [55], a Newton-type penalty method, or any other scheme dedicated to 
such problems. The limit surface can be evaluated at any point using the ideas of Stam [42]. To reduce 
the computational expenses, we apply collision detection on the faceted control mesh instead. There, 
the problem reduces to evaluating the distance between all nine edge-edge pairs and six vertex-face pairs 
belonging to the two triangular elements in question [56]. As a compromise, one could consider doing 
this on any level of subdivision refinement, rather than on the control mesh. 

The closest points on edge-edge pairs are efficiently determined using an algorithm by Sunday [57] 
with edges normalized to unit length. Eriscon's robust algorithm [56] is used for the vertex-triangle pairs. 
Once the closest points have been identified on both contacting triangles, and if they are less than the 
shell thickness apart, a repulsive contact force that is proportional to the surface area of the contacting 
elements is distributed to the six involved vertices using linear interpolation. The precise form of the 
contact law turned out to be irrelevant in many practical applications. An important feature is divergence 
at zero distance to avoid interpenetration when the shell is very thin. 

5 Numerical Examples 

A standard verification obstacle course for small shell deformations in the linear regime can be found 
in Ref. [17]. Although we have verified our implementations using such examples, we are omitting the 
details here, as our focus is on large deformations with geometric nonlinearity, and growth. A few 
tests of nonlinear subdivision shell elements without growth can be found in Refs. [18,45] for different 
constitutive models. In the following, we perform a couple of standard tests to verify our implementation 
of the geometrically nonlinear Koiter shell theory, followed by some examples including various types of 
growth-induced nonlinearity, with increasing sophistication. Unless indicated otherwise, one Gauss point 
is used per element. We denote the subdivision shell finite elements by SD3R in the following for brevity, 
in the spirit of the Abaqus FEA finite element software [58]. 
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5.1 Inflation and Isotropic Growth of a Sphere 



Only few geometrically nonlinear problems are amenable to analytical solution. The inflation of a sphere 
is one of them [59]. For simplicity, we set the bending rigidity D = in this example, i.e., a change 
in energy is assumed purely due to stretching. Since a spherical shell has no boundary, this example is 
perfect for verifying both the pure response to large membrane stresses, and uniform in-plane growth, in 
a single scenario. Consider a growth tensor 

G = diag(l + 0,1 + 0,1) (43) 

with respect to the local tangent basis {ai/|ai|, a 2 /|a 2 |, a 3 }, where g is a positive growth factor. The 
sphere with initial radius R is then trivially expected to grow uniformly according to R/R =: \ = 1 + g. 
On the other hand, in the absence of bifurcations away from the spherical symmetry [60], the pressure 
p needed to inflate a sphere obeying the Koiter energy density W (20) from radius R to R > R is easily 
found by balancing internal and external forces: 

p =7m' (44) 

W = W(R; R) is found using local symmetry on the Green strains 

E ae = \{\ 2 -l)5l (45) 

where A = R/R > 1 is the principal stretch. The energy density of the inflated spherical membrane is 
thus 

W(R;R) = K { -^±(\*-l) 2 , (46) 

yielding the pressure relation 

Yh 

P= W , r(A 3 -A)>0. (47) 

R(l — v) 

The meshes used for this example are shown in Fig. 3. They are constructed by recursive quadrisection 
of the faces of a regular icosahedron, followed by a radial projection of the newly created vertices onto 
the bounding sphere on each level of recursion. The employed recursion depths are 1, 2 and 3, yielding 
triangulated spheres with 80, 320 and 1280 equilateral elements, respectively. The relevant simulation 
parameters are R = 1, h = 10 -3 , Y = 1, v = 0.3. 

In Fig. 4, we plot the growth-expansion and pressure-expansion curves next to each other, demon- 
strating the high accuracy and convergence to the analytical solutions in both cases. Relatively small 
systems are sufficient for high accuracy even at extremely large deformations. 



Figure 3: Icosa-spherical meshes for the inflated and growing spherical shell with 80, 320 and 1280 
triangles (126, 486 and 1926 DOFs). 
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normalized radius A = R/R normalized radius A = R/R 



Figure 4: Change of radius of a spherical shell due to isotropic in-plane growth (left) and uniform pressure 
inflation (right). 



5.2 Pinched Hemispherical Shell 

Next, we turn to two numerical examples without growth, for verification of the linearly elastic, geomet- 
rically nonlinear shell with coupled stretching and bending. The pinched hemisphere is a widely used 
benchmark for "an element's ability to represent inextensional modes" and "rigid body rotations about 
normals of the shell surface" [61]. In its nonlinear regime, it is a test recommended by the National 
Agency for Finite Element Methods and Standards (NAFEMS, 3DNLG-9) [62]. The geometrical setup 
is shown in Fig. 5. A hemispherical shell with an 18° open pole and free boundaries is pinched by four 
equally strong, pairwise opposite diametrical point loads P acting on the equator. The shell radius is 
R = 10, its thickness h = 0.04, and the elastic moduli are determined by Y = 6.825 x 10 7 , v = 0.3. 

In order to minimize the impact originating from the specific choice of boundary constraints, the 
whole hemisphere is simulated without exploiting symmetry, and Green's method is used for the free 
boundaries, i.e., they are unconstrained. In Fig. 6, we plot the displacements of points A and B on the 
limit surface against the applied loads for three mesh resolutions, together with some NAFEMS reference 
results obtained with Abaqus [62,63]. Much more reference data for other nonlinear finite elements can 
be found e.g. in Refs. [62,64,65] and references therein. The employed meshes are obtained by regularly 
discretizing the quarter hemisphere along the angles of inclination and the azimuth, resulting in 128, 512, 
and 2048 triangles per quarter hemisphere, respectively. In Fig. 5, the 16 x 16 mesh with 512 triangles 
per quarter, that is also used in the NAFEMS results, is shown on the right. 

Our load-displacement curves for subdivision surface elements are almost identical to Abaqus S4R. 
High precision is obtained with SD3R at much less DOFs, mainly because subdivision shell elements go 
without rotational variables, unlike all other shown elements. 




Figure 5: Hemispherical shell subject to point loads. Left: Undeformed reference configuration reduced 
to the first quadrant exploiting symmetry. Right: Control mesh of the quasi-static solution at maximum 
load P — 100, without ghost elements. The color encodes the bending energy density. 
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Figure 6: Load-displacement curves for the pinched hemispherical shell. The NAFEMS reference values 
were obtained with Abaqus. 



5.3 Stretched Cylinder with Free Ends 

The next example is another standard loading test, consisting of a cylindrical shell with free boundaries 
that is stretched transversally by two equally strong, opposite diametrical point loads P acting on the 
middle of the cylinder length. This test case has found vast attention in the literature. For an overview, 
see eg. Refs. [64, 65] and references therein. Its peculiar usefulness is due to its ability to examine two 
different response regimes, one after the other. At small loads, the large deformation results from low 
bending stiffness, while at large loads, further deformations require the stiff shell to be stretched primarily. 
The geometrical setup is shown in Fig. 7. The cylinder radius is R = 4.953, its length L = 10.35, its 
thickness h = 0.094, and the elastic moduli are determined by Y = 10.5 x 10 6 , v = 0.3125. Like in 
the previous example, the full shell is simulated neglecting present symmetries, and no constraints are 
applied to the boundaries. In Fig. 8, the resulting load-displacement curves for the points A, B and C 
are compared to Abaqus S4R element data from Ref. [63]. In the bending regime at moderate loads, our 
data from 351 DOFs almost coincides with S4R using 5550 DOFs, again demonstrating the outstanding 
computational efficiency of subdivision shells. In the stretching regime at large loads, SD3R elements are 
slightly less stiff than S4R and elements found in other literature mentioned above. Those displacements, 
however, are quite sensitive to changes in the mesh structure. Other meshes than that shown in Fig. 7 
lead to marginally shifted displacements at large loads. 

Correctly capturing the snap-through transition near P « 2x1 4 has posed a tough challenge to various 
finite shell elements in the past. Some even fail to correctly feature it at moderate mesh resolutions [66-68]. 
With the subdivision shell elements, we have not observed such problems, not even for meshes much 
coarser than mentioned in Fig. 8. 
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Figure 7: Cylindrical shell subject to point loads. Left: Undeformed reference configuration. Middle: 
Control mesh (16x24) of the quasi-static solution at maximum load P = 4xl0 4 , without ghost elements. 
Right: Limit surface of the same configuration, with the stretching energy density on a logarithmic color 
scale. 




concentrated load P x 10 



Figure 8: Load-displacement curves for the stretched cylindrical shell, with the transition near half 
loading. 
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5.4 Anisotropic Growth and Boundary Instabilities 

The purpose of the remaining two numerical examples is to further verify the proper functioning of the 
growing shell elements, and to demonstrate the wide range of possible applications. The first, presented 
in this subsection, addresses an interesting property of plant growth: the development of buckling insta- 
bilities at tissue boundaries, such as the edges of leaves and flower petals. The phenomenon is very similar 
to the permanent deformation left on a plastic sheet that is torn apart. Along the open boundary, wavy 
structures, whose exact nature depends on the shell thickness and growth profile, occur. The out-of-plane 
bending is a result of the tissue's thinness combined with large in-plane strains from compression due 
to in-plane growth (or plastic tension in the tearing case). When a certain critical growth threshold is 
reached, it becomes favorable to bend rather than compress further, and the initial symmetry is broken. 
Consider a shell whose reference configuration is a circular cylinder with radius r(z) = ro, z G [0, L], and 
assume an anisotropic growth tensor 

G = diag(l,l + 0(*),1) (48) 

with respect to the canonical basis of cylindrical coordinates (r, (/?, z). The growth profile g(z) is considered 
monotonically increasing. Following the arguments in Ref. [40], where the Gauss-Bonnet theorem is 
applied, the stability condition for cylindrical symmetry reads 

f z (L)<r , i.e. ^(£)<1- (49) 

We have simulated polynomial growth fields of the form 

g(z) = c(^y, c,p>0, (50) 

for illustration and verification, where the cylindrical symmetry is preserved if c < L/p. A selection of 
fundamental results is shown in Fig. 9. The boundaries are of the free Green type, i.e., the ghost node 
displacements are unconstrained. The essential result of these simulations is that for reasonable angular 
resolutions of the meshed cylinder, the symmetry transition is observed very close to the theoretical value. 
Notice that the linear case p = 1 is equivalent to the excess cone [20,69,70] in the limit ro — » 0, where a 
second phase transition beyond c = L is known to exist. 

For growth profiles with very large gradients (dg/dz)(L) at the boundary, such as 

9(z) =(l + ^) , < l « L, (51) 

where I is a small characteristic length scale, the pattern of the open shell boundary has been reported to 
be self-similar, with an odd integer scaling factor that is mostly 3 and sometimes close to 5 in experiments 
[38,71,72]. This behavior is easily reproducible with our growth implementation. To demonstrate this, 
we have simulated a planar rectangular shell of length L = 1, width W = 4 and thickness h = 10 -4 , using 
a characteristic growth length I = 40ft. The z = edge is clamped using Schweitzer- Cirak constraints, 
while all other edges are unconstrained. The shell is grown using three different mesh resolutions to reveal 
the discretization effects. The employed hierarchical meshes are visualized in Fig. 10. Detailed mesh data 
is listed in Tab. 1, where also the ^-coordinates of the Gauss points closest to the open edge, at which 
the growth tensors are largest, are given for reproducibility. Fig. 11 shows the resulting equilibrated 
deformed configurations after growth. The self-similar nature of the open boundary is discernible even 
on the coarsest mesh, but the finest resolution is needed to get a clear resemblance to experiments [40,71]. 
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growth profile preserved symmetry broken symmetry 




Figure 9: Symmetry breaking for polynomial growth gradients. The reference configuration is cylindrical 
with unit radius and length L = 4. All boundaries are free, and the shell thickness is h = 10 -2 . The 
linear case with 1.5 times the critical growth leads to a ring (ground state: wavenumber k = 2) with 
contact, while the same overcritical growth in the quartic case yields a k = 10 mode at the boundary. 




nodes ghost nodes elements ghost elements shortest element largest quadrature z 

1034 186 1876 376 1/32 0.999457 

4702 590 8804 1188 1/128 0.999887 

38780 4204 73340 8422 1/1024 0.999987 



Table 1: Mesh properties for Figs. 10 and 11. 
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full deformed limit surface control mesh closeup 




Figure 11: Self-similar shell boundary at large growth gradients on different mesh resolutions. All edges 
are free except at z = 0, where the sheet is clamped in all directions. The ghost element layers are not 
shown. In the coarsest configuration (top) already, superimposed waves with scaling factor 3 and 5, as 
predicted in Ref. [38], are observable. The finest resolution (bottom) is detailed enough to manifest about 
three levels of wrinkles, clearly resembling experimentally obtained self-similar wrinkling cascades [40,71]. 



5.5 Confined Growth vs. Crumpling 

We put the growing shell to a final sophisticated test by comparing two processes involving tight self- 
interaction and spatial confinement, in the frictionless, quasi-static, elastic limit. A thin disk with radius 
R is placed inside of a spherical container with the same radius R = R. In the first setup, the container 
is shrunken, and consequently, the disk crumples into a ball of the size of the container, similarly to a 
sheet of paper that one crumples by hand. In the second setup, the opposite happens, i.e., the container 
sustains its size while the disk is subjected to a constant isotropic growth rate, both in plane and in 
thickness. 

Various numerical simulations of crumpled elastic [8,9,73] and elasto-plastic [10] sheets and membranes 
in shrinking spheres have been carried out in recent years. The main finding is that sheets tightly crumpled 
into balls, although consisting mostly of air, develop a very large bulk stiffness resulting from a network 
of ridges and vertices of high magnitudes of mean curvature. A very large portion of the bending energy 
is condensed into this network [5]. A priori, it is not obvious whether a shell growing inside of a fixed 
confinement will exhibit equivalent behavior. With the present thin shell theory, we are able to answer 
this question in the elastic limit. 

The only relevant simulation parameter is the thickness-to-size ratio h/R = 0.01. To obtain equivalent 
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time scales in the two problems, we shrink the sphere in the first case according to R(t) = R/(l + 9(t)), 
where g(t) = Xt is the in-plane growth factor on the growing shell in the second case. Accordingly, to 
achieve equivalent length scales, the growing shell has an increasing thickness h(t) = h(0)(l + g(t)). The 
growth rate A is chosen small enough to allow for a quasi-static simulation in both cases, and damping is 
subcritical. The used mesh consists of 8864 triangles (excluding ghosts, see Fig. 12, first column). 

We have not found any evidence indicating that the two processes are different. During early stages, 
both shells buckle to form a developable cone with a single vertex, where most curvature and bending 
energy is concentrated. Around R/R « 0.53, the apex splits into two vertices, and more vertices subse- 
quently emerge, leading to the same ridge network. In the third column of Fig. 12, the reduced mean 
curvature 

* ■= k ~^R (52) 

is projected onto the unfolded disk, where k\ and are the principal curvatures. The cross correlation 
of the mean curvature ridge patterns is r = 0.89, a very high value when compared to recent similar 
measurements [9]. The fourth column displays the dimensionless rescaled bending energy density 



b 



(53) 



No qualitative or quantitative disparity, going beyond minor local shifts resulting from the finite element 
size and slow but finite dynamics, is observed. While this indicates that the two physical processes are 
in fact equivalent, it also accentuates the large potential of numerical simulations of growing thin shells. 
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Figure 12: Comparison between a crumpled shell (top row) and a growing shell in spherical confinement 
(bottom row). Network vertices are numbered for easy identification. 



6 Conclusions 

We have presented a thin shell finite element approach that is capable of calculating anisotropic in-plane 
growth and large deformations. To this end, the Rodriguez deformation gradient decomposition was used 
to extend the classical Kirchhoff-Love theory. Assuming that the growth tensor is independent of the 
current deformation, the implementation is straightforward and requires only minor modifications, as it 
formally coincides with the implementation of Kirchhoff-Love shells without growth. The presented model 
generalizes the recently popularized target metric approach in that it includes the change of reference 
curvature of the shell. 
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The presented finite elements are of the Loop subdivision surface kind. Subdivision surface interpo- 
lation provides many advantages over traditional C° elements: 

• It meets the continuity requirement imposed by the bending energy and is thus completely locking- 
free and guaranteed to converge. 

• One Gauss point per element is sufficient in many cases. 

• It is very easy to implement and requires only three DOFs per node. No rotational or auxiliary 
variables are required. 

• It can handle arbitrary mesh topologies. 

We have demonstrated the outstanding efficiency of the subdivision shell elements on three standard 
loading examples. Despite their unquestionable superiority to C° elements, subdivision shell elements 
are accompanied by new challenges: 

• The extended local support of subdivision basis functions requires new concepts for boundary 
conditions, adaptive mesh refinement and fracture modeling [74-77]. 

• Contact detection on the limit surface is a nonlinear optimization problem. 

The finite element method is much better suited for strong material anisotropics than many other 
discretization schemes. In this article, anisotropic in-plane growth fields have been built into thin shell 
finite elements without complicating the underlying formalisms. We have illustrated the large range 
of potential applicability of growing finite shell elements to various problems in material science and 
engineering by numerically simulating different growth scenarios. They are also expected to be very well 
suited for the simulation of deformable confining membranes in packing problems [78,79]. A dry friction 
model and plasticity effects will have to be included to simulate real-world time-irreversible contact 
problems. 

This work was financially supported by the ETH Research Grants "Morphogenesis in Constrained 
Spaces" and "Packing of Slender Objects in Deformable Confinements" (ETHIIRA Grants No. TH-06 
07-3 and ETH-03 10-3). 
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